mcmc_intervals2 <- function(data,
transformations = list(),
...,
prob = 0.5,
prob_outer = 0.9,
point_est = c("median", "mean", "none"),
rhat = numeric()) {
# check_ignored_arguments(...)
color_by_rhat <- rlang::has_name(data, "rhat_rating")
no_point_est <- all(data$point_est == "none")
x_lim <- range(c(data$ll, data$hh))
x_range <- diff(x_lim)
x_lim[1] <- x_lim[1] - 0.05 * x_range
x_lim[2] <- x_lim[2] + 0.05 * x_range
# faint vertical line at zero if zero is within x_lim
layer_vertical_line <- if (0 > x_lim[1] && 0 < x_lim[2]) {
vline_0(color = "gray90", size = 0.5)
} else {
geom_ignore()
}
args_outer <- list(
mapping = aes_(x = ~ ll, xend = ~ hh, y = ~ parameter, yend = ~ parameter),
color = get_color("mid")
)
args_inner <- list(
mapping = aes_(x = ~ l, xend = ~ h, y = ~ parameter, yend = ~ parameter),
size = 2,
show.legend = FALSE
)
args_point <- list(
mapping = aes_(x = ~ m, y = ~ parameter),
data = data,
size = 4,
shape = 21
)
if (color_by_rhat) {
args_inner$mapping <- args_inner$mapping %>%
modify_aes_(color = ~ rhat_rating)
args_point$mapping <- args_point$mapping %>%
modify_aes_(color = ~ rhat_rating,
fill = ~ rhat_rating)
} else {
args_inner$color <- get_color("dark")
args_point$color <- get_color("dark_highlight")
args_point$fill <- get_color("light")
}
point_func <- if (no_point_est) geom_ignore else geom_point
layer_outer <- do.call(geom_segment, args_outer)
layer_inner <- do.call(geom_segment, args_inner)
layer_point <- do.call(point_func, args_point)
# Do something or add an invisible layer
if (color_by_rhat) {
scale_color <- scale_color_diagnostic("rhat")
scale_fill <- scale_fill_diagnostic("rhat")
} else {
scale_color <- geom_ignore()
scale_fill <- geom_ignore()
}
ggplot(data) +
layer_vertical_line +
layer_outer +
layer_inner +
layer_point +
scale_color +
scale_fill +
scale_y_discrete(limits = unique(rev(data$parameter))) +
xlim(x_lim) +
bayesplot_theme_get() +
legend_move(ifelse(color_by_rhat, "top", "none")) +
yaxis_text(face = "bold") +
yaxis_title(FALSE) +
yaxis_ticks(size = 1) +
xaxis_title(FALSE)
}n_samp <- NULL
#N <- 2
for (n in 1:length(files)) {
#for (n in 1:N) {
file_name <- str_c("../rda/", files[n])
file_name
load(file_name)
assign(files[n], fit)
n_samp[n] <- list_dat_d$N
}
dat <- tibble(fit = files) %>%
# head(2) %>%
mutate(N = n_samp) %>%
mutate(elpd_list = map(fit, ~ loo(get(.)))) %>%
mutate(looic = map_dbl(elpd_list, ~ .$estimates[3, 1]))
rm(list = ls()[ls() %>% str_detect("\\.rda$")])
DT::datatable(dat)season-habitat-trait_data
files <- list.files("../rda")
#N <- 2
for (n in 1:length(files)) {
#for (n in 1:N) {
file_name <- str_c("../rda/", files[n])
file_name
load(file_name)
tmp_name <- trait6 %>%
select(-sp) %>%
names
tr_name <- c("Int", tmp_name)
title <- str_c(dry, "-", hab, "-", trait_data)
x_name <- names(Xd)
gamma_name <- str_c(
rep(tr_name, length(x_name)),
rep("-", length(tr_name) * length(x_name)),
rep(x_name, each = length(tr_name)))
fig_dat <- mcmc_intervals_data(
fit,
regex_pars = "gamma",
#prob = 0.5,
prob_outer = 0.95,
point_est = "median"
) %>%
mutate(parameter2 = gamma_name) %>%
filter(parameter2 != "Int-Int") %>%
mutate(sig = ifelse(ll * hh < 0, "non-sig", "sig"))
fig_dat
p <- ggplot(fig_dat) +
geom_vline(xintercept = 0) +
geom_segment(data = fig_dat,
aes(x = ll, xend = hh, y = parameter2, yend = parameter2),
col = "darkgrey") +
geom_point(data = fig_dat,
aes(x = m, y = parameter2, col = sig),
size = 4) +
xlab("Coef") +
ylab("") +
ggtitle(title) +
theme_classic()
print(p)
}